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Abstract Gap junctions, also referred to as electrical synapses, are expressed along 
the entire central nervous system and are important in mediating various brain 
rhythms in both normal and pathological states. These connections can form between 
the dendritic trees of individual cells. Many dendrites express membrane channels 
that confer on them a form of sub-threshold resonant dynamics. To obtain insight into 
the modulatory role of gap junctions in tuning networks of resonant dendritic trees, 
we generalise the "sum-over-trips" formalism for calculating the response function 
of a single branching dendrite to a gap junctionally coupled network. Each cell in the 
network is modelled by a soma connected to an arbitrary structure of dendrites with 
resonant membrane. The network is treated as a single extended tree structure with 
dendro-dendritic gap junction coupling. We present the generalised "sum-over-trips" 
rules for constructing the network response function in terms of a set of coefficients 
defined at special branching, somatic and gap-junctional nodes. Applying this frame- 
work to a two-cell network, we construct compact closed form solutions for the net- 
work response function in the Laplace (frequency) domain and study how a preferred 
frequency in each soma depends on the location and strength of the gap junction. 
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1 Introduction 

It has been known since the end of the nineteenth century and mainly from the work 
of Ramon y Cajal [1] that neuronal cells have a distinctive structure, which is dif- 
ferent to that of any other cell type. The most extended parts of many neurons are 
dendrites. Their complex branching formations receive and integrate thousands of in- 
puts from other cells in a network, via both chemical and electrical synapses. The 
voltage-dependent properties of dendrites can be uncovered with the use of sharp mi- 
cropipette electrodes and it has long been recognised that modelling is essential for 
the interpretation of intracellular recordings. In the late 1950s, the theoretical work of 
Wilfrid Rail on cable theory provided a significant insight into the role of dendrites 
in processing synaptic inputs (see the book of Segev et al. [2] for a historical per- 
spective on Rail's work). Recent experimental and theoretical studies at a single cell 
level reinforce the fact that dendritic morphology and membrane properties play an 
important role in dendritic integration and firing patterns [3-5]. Coupling neuronal 
cells in a network adds an extra level of complexity to the generation of dynamic 
patterns. Electrical synapses, also known as gap junctions, are known to be impor- 
tant in mediating various brain rhythms in both normal [6, 7] and pathological [8-10] 
states. They are mechanical and electrically conductive links between adjacent nerve 
cells that are formed at fine gaps between the pre- and post-synaptic cells and per- 
mit direct electrical connections between them. Each gap junction contains numerous 
connexon hemi-channels, which cross the membranes of both cells. With a lumen di- 
ameter of about 1.2 to 2.0 nm, the pore of a gap junction channel is wide enough to 
allow ions and even medium- sized signalling molecules to flow from one cell to the 
next thereby connecting the two cells' cytoplasm. Being first discovered at the giant 
motor synapses of the crayfish in the late 1950s, gap junctions are now known to be 
expressed in the majority of cell types in the brain [11]. Without the need for recep- 
tors to recognise chemical messengers, gap junctions are much faster than chemical 
synapses at relaying signals. 

Earlier theoretical studies demonstrate that although neuronal gap junctions are 
able to synchronise network dynamics, they can also contribute toward the genera- 
tion of many other dynamic patterns including anti-phase, phase-locked and bistable 
rhythms [12]. However, such studies often ignore dendritic morphology and focus 
only on somato- somatic gap junctions. In the case of dendro-dendritic coupling, sim- 
ulations of multi-compartmental models reveal that network dynamics can be tuned 
by the location of the gap junction on the dendritic tree [13, 14]. Here, we develop 
a more mathematical approach using the continuum cable description of a dendritic 
tree (either passive or resonant) that can compactly represent the response of an en- 
tire dendro-dendritic gap junction coupled neural network to injected current using a 
response function. This response function, often referred as a Green's function, de- 
scribes the voltage dynamics along a network structure in response to a delta-Dirac 
pulse applied at a given discrete location. Our work is based on the method for con- 
structing the Green's function of a single branched passive dendritic tree as origi- 
nally proposed by Abbott et al. [15, 16] and generalised by Coombes et al. [17] to 
treat resonant membrane (whereby subthreshold oscillatory behaviour is amplified 
for inputs at preferential frequencies determined by ionic currents such as 4). This 

^ springer 



Journal of Mathematical Neuroscience (2013) 3:15 



Page 3 of 27 




"sum-over-trips" method is built on the path integral formulation and calculates the 
Green's function on an arbitrary dendritic geometry as a convergent infinite series 
solution. 

In Sect. 2, we introduce the network model for gap junction coupled neurons. 
Each neuron in the network comprises of a soma and a dendritic tree. Cellular mem- 
brane dynamics are modelled by an 'LRC (resonant) circuit. In Sect. 3, we focus on 
an example of two unbranched dendritic cells, with no distinguished somatic node, 
with identical and heterogeneous sets of parameters and give the closed form solu- 
tion for network response with a single gap junction. The complete "sum-over-trips" 
rules for the more general case of an arbitrary network geometry are also presented. 
In Sect. 4, we apply the formalism to a more realistic case of two coupled neurons, 
each with a soma and a branching structure. We introduce a method of 'words' to 
construct compact solutions for the Green's function of this network and study how 
a preferred frequency in each soma depends on the location and strength of the gap 
junction. Finally, in Sect. 5, we consider possible extensions of the work in this pa- 
per. 



2 The Model 

We consider a network of cells connected by gap junctions. The neural morphology 
of individual cells includes a branching dendritic structure and a lumped soma (see 
an illustrative example for two cells in Fig. la). We assume that the dendrites are 
not purely passive (i.e. modelled by the 'RC' circuit), but are resonant (i.e. modelled 
by the 'LRC' circuit shown in Fig. lb). Many neurons exhibit resonances whereby 
subthreshold oscillatory behaviour is amplified for inputs at preferential frequencies, 
for example as seen in neurons of rat sensorimotor cortex [18]. In this case, it is 
known that the non-linear ionic current Ih is responsible, and in general it is believed 
that the presence of 4 in dendrites can have a major impact on the integration of 
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subthreshold synaptic activity [19]. From a mathematical perspective, Mauro et al. 
[20] have shown that a linearisation of channel kinetics (for currents such as Ih), 
about rest, may adequately describe the observed resonant dynamics. The resulting 
linear system has a membrane impedance that displays resonant-like behaviour due 
to the additional presence of inductances (which are determined by the choice of 
channel model). This circuit is described by the specific membrane capacitance C, 
the resistance across a unit area of passive membrane R and an inductance L in series 
with a resistance r. The transmembrane voltage Vi (x,t) onan individual branch / of 
each cell is then governed by the following set of equations: 



dVi d^Vi Vi 1 
1 = J). 1 1 

dt dx^ Ti C 



Di^---—[Ii-Iin^,il (1) 



Li ^ = -nli + Vi, 0<x<Ci,t>0. (2) 

The constants Di and Zi can be found in terms of the electrical parameters of the 
cell membrane as Di = ai/(4RajCi) and r/ = Q/?/, where ai is a diameter and 
Raj is the specific cytoplasmic resistivity of branch /. The term li^jj (x,t) models an 
external current applied to this branch. The dendritic structure of each cell is attached 
to an equipotential soma of the diameter as modelled by the 'LRC circuit with the 
parameters Cs, Rs, L^ and r^. Moreover, individual branches of different cells can be 
connected by gap junctions with a coupling parameter Rq]. 

Equations (l)-(2) for each dendritic segment must be accompanied with addi- 
tional equations describing the dynamics of voltage at two ends of a segment. If the 
proximal (x = 0) or distal (x = Ci ) end of a branch is a branching node point the con- 
tinuity of the potential across a node and Kirchoff 's law of conservation of current 
are imposed. For example, boundary conditions for a node indicated in Fig. la take 
the form: 



raj 9-^ 



Vj(Cj,t) = Vn(0,t) = Vk(0,t), (3) 

1 dVn 



ra,n 



(4) 

x=0 



where raj = 4Raj /(naj) is the axial resistance on branch j. If a branch terminates 
at X = jCi we either have a no-flux (a closed-end) boundary condition 



dx 



= 0, (5) 



or a zero value (an open-end) boundary condition 

Vi(Ci,t) = 0. (6) 

A lumped soma can be treated as a special node point with the somatic membrane 
voltage Vs(t) and the following set of equations, which imposes special boundary 
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conditions on the proximal ends of branches connected to the soma: 



V,{t) = Vj{0, t), 



dV^ 
df 



1 dVi 



dx 



-h 



x=0 



(7) 
(8) 

(9) 



where the sum in Eq. (8) is over all branches connected to the soma. If the branches 
of two cells are coupled by a gap junction, the location of this coupling can be treated 
as a special node point on an extended branching structure. This gap-junctional (GJ) 
node requires the following set of boundary conditions (given here with an assump- 
tion that it is placed at x = 0): 



Vm-(0,t) = VmAO,t), 



(10) 



and 



1 /dV^ 



dx 



1 /aK 



+ 



x=0 



dx 



^ a,n \ 



x=0 9X 



x=0 



x=0 



■gGj{Vn,-(0,t)-V^-(0,t)), (11) 
■gGj{Vn-(0,t)-V^-(0,t)), (12) 



where ^gj = V^gj is the conductance of the gap junction and m~ and (n~ and 
n^) are two segments of branch m (branch n) connected at the gap junction (see 
Fig. la). The expressions in (10) reflect continuity of the potential across individual 
branches m and n, and Eqs. (11)-(12) enforce conservation of current. 

A whole network model can be viewed as an extended tree structure with each 
individual node belonging to one of the following categories: a terminal, a regular 
branching node, a somatic node or the GJ node. The voltage dynamics along the net- 
work structure are described by linear equations and, therefore, the model's behaviour 
can be studied by constructing the network response function known as the Green's 
function, Gij(x, y,t). This function describes the voltage response at the location x 
on branch / in response to a delta-Dirac pulse applied to the location y on branch 
j at time ^ = 0 (branches / and j can belong either to the same cell or to the two 
different cells). Knowing the Green's function for the whole structure, it is easy to 
compute the voltage dynamics along the whole network for any form of an external 
input /inj,j(-^, 0 applied to branch j as 

Vi{x,t) = y] [ ' dyGik(x,y,t)Vk{y,0) 

ds dyGij(x,y,t -s)Iinjj(y,s), 
Jo Jo 



+ 



(13) 
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Fig. 2 A network of two 
identical cells, each consists of 
an infinite dendritic cable, 
coupled by a gap junction 



-^inj,m 



>^ m 



i) y X 



x = 0 



D 



D 



where Vk (x , 0) describes the initial conditions on branch k and the sum is over all 
branches of the tree. Multiple external stimuli can be tackled by simply adding new 
terms with additional inputs Iinjj(x, t) to Eq. (13). 



3 The Green's Function on a Network 

Earlier work of Coombes et al. [17] demonstrated that the Green's function for a sin- 
gle cell with resonant membrane can be constructed by generalising the "sum-over- 
trips" framework of Abbott et al. [15, 16] for passive dendrites. Here, we demonstrate 
how this framework can be extended to a network level starting with the simple case 
of two identical cells. 

3.1 Two Simplified Identical Cells 

We consider the case of two identical cells coupled by a gap junction. Each cell 
consists of a single resonant dendrite of infinite length (see Fig. 2). A gap junc- 
tion controlled by the parameter Rqj and located at x = 0 divides two dendrites into 
four semi-infinite segments: m~,m^,n~, and n~^. We assume that an external input 
^inj,m- (x,t) = 8(x — y)8(t) is applied to segment m~ . The Green's function on each 
segment must satisfy the set of Eqs. (l)-(2) with the boundary conditions at the gap 
junction given by Eqs. (10)-(12). Introducing the Laplace transform with spectral 
parameter co 



/'OO 

4/(0] = /M= / e-"7(0 
Jo 



dt, 



and assuming zero initial data, we can solve this model in the frequency domain: 



G^-(x,y,co) = 
G^+(x,y,a)) = 



■ Q-y{co)\x-y\ 
2Dy(co) 

PG](co))- 



PG](o^)- 



-y{co)\x-\-y\ 



2Dy(co) 

-y{a))\x-\-y\ ' 



2Dy(co) 



Gn-(x,y,a)) = Gn+(x,y,a)) = 



-y(aj)\x-\-y\ " 

2Dy(co) 



(14) 
(15) 
(16) 
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Fig. 3 All trips (including the 
shortest trip for G^- {x,y,co) 
shown by a dashed arrow) and 
the corresponding coefficients 
Atrip ((^) at the GJ for two 
identical cells 



-PGJ 



x = 0 



D 



1 -PGJ 
- PGJ 



X 



D 



where 

in 11 

+ ^+T7— 7T , (17) 



9 1 
y^{a)) = — 



r C(r-\-coL) 



and 



1 

^gjM = ;777-^^^ TTT' z(co) = y(co)/ra. (18) 

2{z(co)Rgj + 1) 

Solutions (14)-(16) are obtained using the "sum-over-trips" method where G(x,y,co) 
on each segment can be found as J^^j-ips ^trip(<^)^oo(Arip, <^), and 

^ Q-y(a))\x\ 

G^(x,co) = ------ (19) 

2Dy((jo) 

is the Laplace transform of the Green's function Goo(x,t) for an infinite resonant 
cable. £trip is the length of a path that starts at point x on one of the segments and 
ends at point y on segment m~ . The trip coefficients Atri^ico) which ensure that the 
boundary conditions at the gap junction hold are chosen according to the following 
rules (see Fig. 3): 

• Atnp(co) = —pG]((^) if the trip reflects along on the gap junction back onto the 
same dendrite. 

• Atrip ((w) = 1 — PG]((^) if the trip passes through the gap junction along the same 
dendrite. 

• Atnp(co) = PG]((^) if the trip passes through the gap junction from one cell to 
another cell. 

Performing the numerical inverse Laplace transform (£~^) of Eqs. (14)-(16), we 
obtain the Green's function in the time domain for each segment. These Green's func- 
tions are plotted in Figs. 4a-c. For any arbitrary form of external input Im],m-(^^0 = 
8(x — y)I(t), the voltage response on each segment can be found by taking a con- 
volution of the corresponding Green's function with this stimulus. Using the Laplace 
representation of the Green's function on each segment given by Eqs. (14)-(16) this 
can be computed as 

Vk(x, t) = [Gk(x, y, co)T(co)], k e {m-,m^, n-,n^}, (20) 

where I(co) = £[I(t)]. In Figs. 4d-f, we plot the voltage profiles on each segment in 
response to a rectangular pulse I(t) = r]{)G{t)0{TR — t) of strength ?7o and duration 
xr, where G(t) is the Heaviside step function. 
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0 2 4 6 8 10 12 0 2 4 6 8 10 12 ^0 2 4 6 8 10 12 

Time (ms) Time (ms) Time (ms) 




0 10 20 30 40 50 0 10 20 30 40 50 0 10 20 30 40 50 
Time (ms) Time (ms) Time (ms) 



Fig. 4 a-c The Green's functions G^- (x, y, t), (x, y, t), and G^± (x, y, t) for a model in Fig. 2 
when X = 10 |um and y = 100 |um. Parameters: a = 2 |um, D = 50000 |um^/ms, t = 2 ms, C = 1 |uF/cm^, 
= 100 Q cm, r = 100 Q cm^, L = 5 Hcm^, Rqj = 100 MQ. d-f Voltage profiles on each segment in 
response to a rectangular pulse of strength r]Q = 2nA and duration t^^ = 5 ms applied to segment m~ . 
Note that different y-axis limits are used in c and f 

A response of the network model is characterised by the Green's function and 
can be studied by introducing a power function Vk(x,y,(jo) defined SisVk(x,y,(jo) = 
\Gk(x, y, co)]'^. Resonant dynamics of the model for a given pair of locations (x, y) 
are directly linked with a value f2o at which the function Vk(x,y,(jo) has its max- 
imum. In Figs. 5a-c, we plot the voltage profiles on each segment in response to 
a chirp stimulus /chirp (0 = ^chirp sin(a;chirp^^). These figures clearly demonstrate 
resonant behaviour of the system maximising the voltage responses for particu- 
lar frequencies. In Fig. 5d, we plot the normalised power functions (x, y,co) = 
Vk(x,y, co)/m2iXa)[Vk(x, y,co)] at the same locations. These power functions have 
their maximum at the value f2o = 0.4271, the same for each segment and, therefore, 
the resonances (indicated by arrows) in Figs. 5a-c occur at the same time. We can 
also notice that the function (x,y,(jo) decays to zero more rapidly than the func- 
tions Vj^-ix, y, co) and Vj^+ix, y,co). This explains the rapid reduction of voltage 
amplitude straight after the resonance in Fig. 5c. 

Typical values of a unitary gap junction conductance are 10-550 pS [11] giving 
= 10-5500 pS (or Rqj = 180-10^ M^2) for 1-10 gap junction channels per elec- 
trical connection, although these estimates may be conservative and conductances 
from a larger range could be considered, as for example ^gj = 10-240000 pS cor- 
responding to Rq] = 4-10^ M^2 in [14]. To demonstrate how the resistance of the 
gap junction affects the response function in a model of two identical cells, we plot 
Fig. 6 for Rq] = 100 M^2 (black curves, the case shown in Figs. 4a-c), Rq] = 1 M^2 
(red curves) and Rq] = 1000 MQ (green curves). Low gap-junctional resistance sig- 
nificantly increases the amplitude of the Green's function for Cell n and reduces the 
amplitude of the Green's function for Cell m. Increasing the resistance reduces the 
response in Cell n and slightly increases the response in Cell m . In this model of two 
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0 50 100 150 200 250 300 0 2 4 6 8 10 

Time (ms) ^ 



Fig. 5 a-c Voltage profiles on each segment at x = 10 jam in response to a stimulus /chirp (0 applied at 
y = 100 |im on segment m~ . Cells' parameters as in Fig. 4, c^chirp = 0.003, A^hirp = 1 nA. d Normalised 
power functions 7^^_ (x, y, co) (black curve), 7^^^ (x, y, co) (red curve), (x, y, co) (green curve) 




02468 10 12 02468 10 12 02468 10 12 



Time (ms) Time (ms) Time (ms) 

Fig. 6 a-c The Green's functions G^- (x, y, t), (x, y, t), and G^± (x, y, t) for a model in Fig. 2 
when :v; = 10 |um and y = 100 |um. Parameters: a = 2 \xm, D = 50000 |um^/ms, t = 2 ms, C = 1 |iF/cm^, 
Ra = 100 Qcm, r = 100 ^^cm^, L = 5Hcm2, Rqj = 100 (black curves, as in Figs. 4a-c), 
Rqj = 1 MQ (red curves), Rq] = 1000 MQ (green curves). Note that different y-axis limits are used 
in c 

identical cells, the change of the resistance of the gap junction does not affect the 
resonant frequency which is the same for each segment. 

Although solutions (14)-(16) are found for the case of the resonant membrane, an 
'LRC circuit can be naturally turned into a 'RC circuit by using the limit r ^ oo 
which gives y^ico) = (1/r -\- <jo)/D. In the case of purely passive membrane, it is 
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Time (ms) Time (ms) Time (ms) 



Fig. 7 a-c The Green's functions G^- (x, y, t), (x, y, t), and G^± (x, y, t) for a model in Fig. 2 
with passive membrane when x = 10 |um and y = 100 |um. Parameters: a = 2 |um, D = 50000 |im^/ms, 
r = 2 ms, C = 1 |uF/cm^, Ra = 100 Qcm, Rqj = 100 MQ. d-f Voltage profiles on each segment in 
response to a rectangular pulse of strength = 2 nA and duration tj^ = 5 ms applied to segment m~ . 
Note that different y-axis limits are used in c and f 



possible to make extra progress and find analytical forms of the solutions in the time 
domain (see Appendix A). In Figs. 7 a-c, we plot the Green's functions for the model 
in Fig. 2 with passive (instead of resonant) membrane. Voltage responses on each 
segment in response to a rectangular pulse are shown in Figs. 7d-f. 

3.2 Two Simplified Non-identical Cells 

Here, we consider a model in Fig. 2 with the assumption that the cells are non- 
identical. Then using the Laplace transform and solving the model (with zero initial 
data) in the frequency domain, we obtain 



Gm-(x,y,co) 
G^+(x,y,co) 



1 

Cjfi 



Cff 



' Q-ym(co)\x-y\ 
(1 - PGJ,n(^)) 



PG],n(0^)- 



-ymico)\x-\-y\ 



-ymM\x-\-y \ ■ 



2Z)mKm(<w) _ 



Gn-(x,y,co) = Gn+(x,y,a)) = — 



PGJ,m(co)- 



-\yn{co)x-\-ymico)y\ " 
2Dtnym(co) 



(21) 
(22) 
(23) 



where the parameters Ymioj) and Ynicjo) are defined in terms of cells' individual prop- 
erties as 



1 

1 



■ 1 

1 



- CO + 



1 



1 

C„(r„ + wL„) 



(24) 
(25) 
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Fig. 8 All trips (including the 
shortest trip for G^- {x,y,co) 
shown by a dashed arrow) and 
the corresponding coefficients 
Atrip ((^) at the GJ for two 
non-identical cells 



-PGJ,n 



d 



x = 0 



Cell m 



1 -PGJ,n 
- PG3,rj 



D 



X 



Celln 



Solutions (21)-(23) show that the trip coefficients Atrip (o^) depend on either pGJ,m 
or PGin((^) (see Fig. 8), which have the forms 

PGlmico) = -— , TT' = KmM/^a,m, (26) 

Zm (<^) + Zn (m) + 2RG]Zm ((^)Zn (o)) 

PGln(0)) = -— , , X , Zn(c^) = yn{co)/ra,n- (27) 

Zm (<W) + Zn {co) + 2RGJZm ((^)Zn (co) 



In Figs. 9a-c, we demonstrate how individual variations in cell parameters affect 
the voltage response in the system. For each set of the parameters, we plot the Green's 
functions G^- (x,y,t), (x,y,t), and G^± (x,y,t) obtained by taking the numer- 
ical inverse Laplace transform of (21)-(23). Black curves show the profiles for two 
identical cells. Dashed red curves are the Green's functions for a case when Ln is 
changed from 5 Hcm^ to 25 Hcm^. This change affects the response in Cell n, but 
not in Cell m. Blue curves are plotted for a case when is changed from 5 Hcm^ 
to 1 Hcm^. It has noticeable effect on both cells. Finally, green curves are plotted 
for the case am = 1 jim instead of the original diameter =2 |im showing changes 
in profiles in both cells. As the stimulus in these examples is applied to Cell m, any 
variations in the parameters of this cell have an immediate effect on the responses in 
Cell n. In contrast. Cell m seems to be mostly robust to variations in parameters in 
Cell n. Resonant properties of the cells' responses can be studied by plotting the nor- 
malised power functions V^(x, y, co) for each of the parameter sets (see Figs. 9d-f). 
The heterogeneity of the cells' parameters leads to appearances of different values of 
Qq (the maximum of the power function) for each cell. We can also notice that the 
power functions for Cell n are more localised around their peaks (Fig. 9f) in compar- 
ison to the power functions for Cell m (Figs. 9d, e) as it has been earlier observed in 
the case of two identical cells. 



3.3 An Arbitrary Network Geometry 



Here, we consider a network of spatially-extended cells of arbitrary geometries. This 
network can be treated as a single extended tree structure which consists of individual 
branches (indexed by a finite non-repeating sequence {1, 2, ^, j, .. .}) 
and three types of nodes: a regular branching node, a somatic node and a GJ 
node (see Fig. 10). The Green's function Gij(x,y,t) for the whole structure can 
be found by obtaining Gij(x,y,co) in the Laplace domain and then performing 
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Fig. 9 The Green's functions G^- (x, y, t), (x, y, t), and G^± (x, y, t) (a-c) and normalised power 
functions 7^^_ {x,y,(jo), 7^^^ (x,y,(jo), and (x,y,(jo) (d-f) for a model in Fig. 2 when :sc = 10 |um 
and y = 100 jam. curves: two identical cells with the parameters as in Fig. 4. Dashed red curves: 

as in an identical case except L„ = 25 Hcm^. 5/^^ curves: as in an identical case except Lm = 1 Hcm^. 
Green curves: as in an identical case except flm = 1 l^m 



Fig. 10 A network of two cells 
as an extended tree structure 
with different types of nodes 



Somatic node 



Somatic node 




Terminal 



Brandling node 



£^~^[Gij(x, y, co)]. We consider a general case when each branch of the network 
can have different biophysical parameters and is characterised by the function ykico) 
defined as 



1 



1 1 

\-co-{ 

n Ck(rk + o)Lk) 



(28) 



where k labels an arbitrary branch of the network. Using the "sum-over- trips" for- 
malism Gij(x, y, co) can be constructed as an infinite series expansion 



Gij(x,y,co) 



1 



DjYjico) 



^ Atrip (6t>)Hoo(AripO', j,X,y,C0)), (29) 



trips 



where Hoo(x) = e '^'/^ and >Ctrip(/, j,x,y,co) is the length of a path along the net- 
work structure that starts at the point yi (co)x on branch / and ends at the point yj (co)y 
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on branch j . Note that the length of each branch of the network needs to be scaled 
by yk((^) before >Ctrip is calculated for (29). It is also worth mentioning here that if 
all branches of a network have the same biophysical parameters, i.e. yk(co) = y(co), 
the function Hoo (Ctrip (co)) / (Dy (co)) = Goo(Arip, oj) defined by (19). The trip coef- 
ficients Atrip ((^) in (29) are chosen according to the following set of rules: 

• Initiate Atnp(co) = 1. 

Branching node 

• For any branching node at which the trip passes from branch / to a different 
branch k, Atrip (a;) is multiplied by a factor 2pk((jo). 

• For any branching node at which the trip approaches a node and reflects off this 
node back along the same branch k, Atrip ((w) is multiplied by a factor 2pk(co) — 1 . 

Here, the frequency dependent parameter pk (co) is defined as 

Pk (co) = — , Zk (co) = , (30) 

2^n^n(co) ra,k 

where the sum is over all branches connected to the node. 
Terminal 

• For every terminal which always reflects any trip. Atrip is multiplied by +1 for the 
closed-end boundary condition or by — 1 for the open-end boundary condition. 

Somatic node 

• For the somatic node at which the trip passes through the soma from branch / to 
a different branch k, Atrip (a;) is multiplied by a factor 2ps^k(co). 

• For the somatic node at which the trip approaches the soma and reflects off 
the soma back along the same branch k, Atnp(co) is multiplied by a factor 

^Ps,k(^) - 1- 

Here, the frequency dependent parameter Ps,k(co) is defined as 

Zk(co) 1 1 

Ps k(co) = , ys(co) = CsCO -\ \ , (31) 

2_^nZn(co) + ys(co) Rs rs + LsCO 

where the sum is over all branches connected to the soma. 
GJ node 

• For the GJ node at which the trip passes through the gap junction from branch n 
to branch m. Atrip ((w) is multiplied by a factor pG],m (co). For the GJ node at which 
the trip passes through the gap junction from branch m to branch n, Atnp(co) is 
multiplied by a factor PGJ,n(co). 

• For the GJ node at which the trip approaches the gap junction, passes it and 
then continues along the same branch m, Atnp(co) is multiplied by a factor 
1 — PGin(co). For the GJ node at which the trip approaches the gap junction, 
passes it and then continues along the same branch n. Atrip (a;) is multiplied by a 
factor 1 - pGJ,m(^)' 

^ springer 



Page 14 of 27 



Y. Timofeeva et al. 



Fig. 11 A two-cell model Cell 1 Cell 2 




• For the GJ node at which the trip approaches the gap junction and reflects off 
the gap junction back along the same branch m, Atrip ((w) is multipHed by a factor 
—PGin(o^)- For the GJ node at which the trip approaches the gap junction and 
reflects off the gap junction back along the same branch n, Atrip (a;) is multiplied 
by a factor -pGimif^)- 

Here, parameters PG],m((^) and PG],n((^) are defined by Eqs. (26) and (27). 

We refer the reader to Coombes et al. [17] for a proof of rules for branching and 
somatic nodes. In Appendix B, we prove that the rules for generating the trip coeffi- 
cients at the GJ node satisfy the gap-junctional boundary conditions. 

4 Application: Two-Cell Network 

Here, we demonstrate how the "sum-over- trips" formalism can be applied to a two- 
cell network for obtaining insight into network response. As an example, we consider 
a model of two identical cells, each of which consists of a soma and N attached semi- 
infinite dendrites as shown in Fig. 11. The cells are coupled by a dendro-dendritic 
gap junction located at some distance £gj away from their cell bodies. We assume 
that this network receives an input at the location yo. To study the dynamics of this 
network, we use the "sum-over-trips" framework and construct the Green's functions 
Gi(x, yo, co) and G2(x, yo, co) for Cell 1 and Cell 2, respectively. 

4. 1 Method of Words for Compact Solutions 

Here, we introduce a method which allows us to construct compact solution forms for 
the Green's functions of this two-cell network. We describe this method in detail by 
constructing the Green's function G2(xo, yo, <^) for Cell 2 when xo is placed between 
the soma and the gap-junction as shown in Fig. 12. Introducing points from 1 to 4 on 
this network, we associate letters with different directions as follows: 

• From xo ^ 1 or from 2—>\ \ letter A. 

• From Xo ^ 2 or from 1^2: letter B. 

• From 3^4: letter W. 
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Fig. 12 A two-cell model with 
associated letters 



Cell 1 



Cell 2 




• From 4 

• From 3 



3: letter Y. 
yo: letter Z. 



Then the shortest trip which is a trip from ^ 2 ^ 3 ^ is associated with the 
(ordered) word BZ consisting of one syllable. This and any trip associated with a 
word that starts with the letter B, i.e. from ^ 2, and ends with the letter Z, i.e. 
from 3 ^ yo, will belong to class 1. Any trip associated with a word which starts 
from xo ^ I ^ 2 and ends with 3 ^ yo, will belong to class 2. The shortest trip 
in this class is associated with the word ABZ, consisting of two syllables, AB and 
BZ. We introduce the following table that associates individual coefficients in the 
"sum-over-trips" framework with the syllables: 





A 


B 


w 


Y 


z 


A 


0 


2ps(co) - 1 


0 


0 


0 


B 


-PG](C0) 


0 




0 




W 


0 


0 


0 


2ps(co) - 1 


0 


Y 




0 




0 




Z 


0 


0 


0 


0 


0 



As the cells are identical in this network, the parameter /?gj(<^) is defined by (18) and 

y(co)/ra 



Ps(0)) ■ 



Ny(a))/ra + ys(o))' 



(32) 



where N is 3. number of dendrites attached to each soma and Ysico) is given in (31). 
Then, using the table, it is easy to conclude that for example, the word BZ is as- 
sociated with the coefficient /?gj(<^) and the word ABZ, consisting of the sylla- 
bles AB and BZ, is associated with the coefficient (2ps(co) — l)pG](co)- We also 
notice from the table that different coefficients are associated with the syllables 
BZ and YZ and, therefore, we need to introduce two additional classes. Class 3 
will include the trips with the main skeleton xo^2^3^4^3^yo and 
the associated word BWYZ. Class 4 will include the trips with the main skeleton 
xo^l^2^3^4^3^3;o and the associated word ABWYZ. Combining the 
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skeleton structures (the shortest words) of the four classes, we have 



BZ + ABZ + BWYZ + AB WYZ = (1+A)BZ +(1+A)BWYZ. 



(33) 



class 1 and class 2 class 3 and class 4 



Any new word in each class can be formed by adding a combination of syllables 
AB and WY into the structure (33). Introducing such additions of n combinations of 

syllables consisting of k syllables AB and (n — k) syllables WY by [ ] (both 

syllables AB and WY can take any position in this sequence of n syllables), class 3 
and class 4 can be generalised as 



Similarly, we can generalise class 1 and class 2. However, to ensure that the words 
belong to class 1 and class 2, the syllable AB must be at the end of each word in 
combinations. This can be written as 



(1+A)B[ ]WYZ. 

^ ' 



(34) 



class 3 and class 4 



(1+A)B[ ]'Z, 



(35) 



class 1 and class 2 



where 



[ ]' = [! + [ ]AB]. 



(36) 



Using combinatorics, we can write 





which leads to 




n-l 




■)) 



n-k-1 





^ Springer 



Journal of Mathematical Neuroscience (2013) 3:15 



Page 17 of 27 



+ (^^{2ps(co) - iy{-PGj(co)y~\ (38) 

Considering the possible trips in class 3 given by 

B[ ]WYZ (39) 

and substituting the expression for [ ] found in (38), we find 

PGJ(0))(^^ {2ps((^) - iy{-pGi(co)y~\-pGj((^)){2ps((^) - l)(l - PGiico)) 

n-l y X 

- (-m(^))l] Q(2/^s(^) - iy{-PGj(co)y~^ 

X {-pG](co)){2ps(co) - l)(l - PG](C0)) 

+ {-PGj(c^))(^^{2ps(c^) - iy{-PGj(c^)y~^ 

X PGj(co){2ps(co) - l)(l - PGJ(C0)) 

= [2^{2ps(w) - iy{-pGj(o^)y]pGj((^){2ps(o^) - 1)(1 - PGj(co)). (40) 
Similarly, the possible trips in class 4 given by 

AB[ ]WYZ (41) 

generate the trip coefficients 

{2ps(w) - l)[2%2ps(w) - iy{-pGj(co)y]pGj(co){2ps(w) - l)(l - m(^)). (42) 

Expressions (40) and (42) for coefficients in the trips belonging to class 3 and class 4 
can now be used in the "sum-over- trips" expansion (29) to obtain 

00 

Y^[-2pG](co){2ps(co) - l)YpG](co){2ps(co) - l)(l - PGj((^)) 
X [Goo{yo -xo + 2(n + 1)Cgj,(jo) 

+ {2ps(a)) - l)Goo{yo + 2(n + 1)/:gj, c^)\ (43) 

where Goq{x,o)) — Hoo(y ((jo)x) / (Dy (co)) defined by (19). Similarly, we can show 
that the possible trips in class 1 and class 2 given by 

(l+A)B[l + [ ]AB]Z (44) 
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generate the terms 

PG](co)Goo(yo-XO,Oj) + pG](co){2ps(co) l) Goo (^0 + -^0, 

oo 

+ Y^[-2pG](co){2ps(co) - l)Y PG](co){-PG](co)){2ps(co) l) 

n=0 

X [Goo{yo + 2(n + 1)£gj, 

+ {2p,(co) - l)Goo{yo + 2(n + 1)£gj, ^)]. (45) 
Combining together (43) and (45), we obtain 

<^2(-^0, yo, ^^) = m(<^)<^oo(yO - ^0, ^) + J9Gj(<w)(2j9s(6t>) - l)Goo(yo + -^0, ^) 

oo 

+ ^[-2m(^)(2;^s(^) - i)]Vgj(^)(2/^s(^) - 1) 

X {(-m(<^))[<^oo(yo - -^0 + 2(^ + 1)>Cgj, 

+ i^p^icD) - l)Goo(yo + -^0 + 2(/i + 1)/:gj, ^^)] 

+ (l - m(<^)) [<^oo(yo - ^0 + 2(n + 1)>Cgj, 0^ 

+ (2/?s(^^) - l)Goo(yo +-^0 + 2(/t + 1)/:gj,^)]}. (46) 

Terms in {• • • } in (46) represent multiple trips in each of four classes and since there 
is a match in the length of trips among different classes, Eq. (46) can be simplified as 

G2(x^,y^,(D) = pG](co)[Goo(yo - ^0, oj) + {2ps(co) - l)Goo(yo +-^0,^)] 

oo 

+ ^2^(-/7Gj(^)(2/^s(^) - l)y^\2pGj(co) - 1) 

n=0 

X [Goo(yo - -^0 + 2(n + 1)£gj, ^) 

+ {2ps(co) - l)Goo(yo + -^0 + 2(n + 1)Cgj, co)]. (47) 

Using this method of 'words', we can construct compact solution forms for the 
Green's function for each of these two cells for any combinations of input, y, and 
output, X, locations. For example, placing xq in Cell 1 between its soma and the 
gap-junction we obtain 

Gi(xo, yo,oj) = (l - PG](co))[Goo(yo -xo,co) + {2ps(co) - l)Goo(yo + -^0,^^)] 

oo 
n=0 

X [Goo(>'o -X0 + 2(n + 1)£gj, 

+ (2/?sM - l)Goo{yo+xo + 2{n + 1)£gj, «)]• (48) 
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0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500 

Fig. 13 Values of ^2q at the soma of each cell as a function of Cqj when yQ = jCqj + 10 |um and 
Rqj = 100 MQ (red circles), Rqj = 1000 MQ (black crosses). Passive somas with parameters: diame- 
ter as = 25 |um, Cs = 1 |uF/cm^, = 2000 Q cm^. Resonant dendrites (N = 4) with parameters given in 
Fig. 4 

4.2 Network Dynamics 

To study the role of a gap-junction in this two-cell network model, we focus on the 
Green's functions at the somas of these two cells in response to a stimulus at the loca- 
tion yo. Using Eqs. (47) and (48), we obtain the following somatic response functions: 

Gi(0,yo,co) = 2ps(co){l - PG](co))Goo(yo, co) 

oo 

+ ^2^(-m(^)(2/^s(^) - l)y^\l-2pGAco)) 

X 2ps(co)Goo{yo + 2(n + l)jCGh co), (49) 

and 

G2(0, yo, co) = 2ps(co)pG](co)Goo(yo, co) 

oo 

+ ^2^(-/7Gj(a;)(2/7s(a;) - l)y^\2pGj(co) - l) 

n=0 

X 2ps(co)Goo{yo + 2(n + 1)Cgj, co). (50) 

Resonant properties of each cell are analysed by studying a preferred frequency f2o 
for each cell. This is defined as the frequency at which the corresponding power 
function, Viio)) = |Gi(0, yo, a))\^ for Cell 1 and V2(o)) = 1^2(0, yo, a))\^ for Cell 2, 
reaches its maximum. This means that f2o for each soma is simply a solution of one 
of the corresponding equations, dVi((jo)/dco = 0 and dV2{oo>) /dco = 0. 

In Fig. 13, we plot how this preferred frequency ^o varies as a function of the 
distance of the gap-junction away from each soma. These plots are obtained for the 
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0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500 



Fig. 14 Values of ^2q at the soma of each cell as a function of jCqj when yQ = Cqj + 10 |um and 
Rqj = 100 (red circles), Rqj = 1000 MQ (black crosses). Resonant somas with parameters: diam- 
eter fls = 25 |um, Cs = I |uF/cm^, /?§ = 2000 Qcm?', rs = I Qcw?', Lg =0.1 Hcm^. Passive dendrites 
(N = 4) with parameters given in Fig. 4 




0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500 

Fig. 15 Values of at the soma of each cell as a function of jCqj when yQ = Cqj + 10 |um and 
Rqj = 100 M^2 (red circles), Rqj = 1000 MQ (black crosses). Resonant somas with parameters: diame- 
ter as = 25 [im, Cs = I |iF/cm^, /?§ = 2000 Q cm^, rs = I Q cm^, Ls =0.1 Hcm^. Resonant dendrites 
(N = 4) with parameters given in Fig. 4 



case of a passive soma and resonant dendrites. In Fig. 14 we demonstrate how f2o is 
affected if a resonant soma and passive dendrites are considered. Finally, in Fig. 15, 
we plot f2o(CG]) when each soma and dendritic branch is modelled with resonant 
membrane. All these three figures clearly demonstrate that the somatic response in 
each cell strongly depends on the location of the gap-junction and it is tuned to be 
maximised for different frequencies. This can be shown by applying a chirp stimulus 
at the location yo and plotting the somatic voltage for each cell (see Fig. 16). Reso- 
nances in each cell occur at different times as predicted by Fig. 15. All these figures 
are constructed for truncated series solutions (49), (50) when the index n increases up 
to 20, although it is possible to show that the solutions rapidly converge for a much 
smaller n such as /t = 10. 
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Time (ms) jime (ms) 




0 100 200 300 400 500 0 100 200 300 400 500 

Time (ms) Time (ms) 

Fig. 16 Voltage profiles in the somas of cells in response to a stimulus /chirp (0 applied at the location 
yo = ^GJ + 10 l^^- Cells' parameters as in Fig. 15, Rqj = 100 MQ, (Wchirp = 0.003, A^hirp = 1 nA. Black 
curves: Cqj = 50 |um, green curves: Cqj = 500 |um 

5 Discussion 

In this paper, we have generaUsed the "sum-over-trips" formaUsm for single dendritic 
trees to cover networks of gap-junction coupled resonant neurons. With the use of 
ideas from combinatorics, we have also introduced a so-called method of 'words' 
that allows for a compact representation of the Green's function network response 
formulas. This has allowed us to determine that the position of a dendro-dendritic 
gap junction can be used to tune the preferred frequency at the cell body. Moreover 
we have been able to generate mathematical formula for this dependence without 
recourse to direct numerical simulations of the physical model. One clear prediction is 
that the preferred frequency increases with distance of the gap junction from the soma 
in a model with passive soma and resonant dendrites. In contrast for a system with a 
resonant soma and passive or resonant dendrite, the preferred frequency decreases as 
the gap junction is placed further away from the cell body. 

There are a number of natural extensions of the work in this paper. One is an 
application to more realistic network geometries or more than just two neurons, as 
may be found in retinal networks. Here, it would also be interesting to exploit any 
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network symmetries (either arising from the identical nature of the cells, their shapes, 
or the topology of their coupling) to allow for the compact representation of network 
response (and further utilising the method of 'words'). Another is to incorporate a 
model of an active soma whilst preserving some measure of analytical tractability. 
Schwemmer and Lewis [21] have recently achieved this for a single unbranched cable 
model by coupling it to an integrate-and-fire soma model. The merger of our approach 
with theirs may pave the way for understanding spiking networks of gap junction 
coupled dendritic trees. Moreover, by using the techniques developed by them in 
[22] (using weakly coupled oscillator theory) we may further shed light on the role 
of dendro-dendritic coupling in contributing to the robustness of phase-locking in 
oscillatory networks. 
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Appendix A: Two Simplified Identical Cells with Passive Membrane 

Equations (14)-(16) with y^(co) = (l/r -\- co)/ D provide the solutions of a model in 
Fig. 2 with passive membrane. We introduce the function 



F(x, CO, q) ■■ 



1 



-yM\x\ 



y{co)-\-q 2Dy((jo) 

and its inverse Laplace transform 

F(x,t,q) = ^~^[F(x,co,q)] 

= lel-l^e^^'^-i/^)^erfcf^yD + ^ 

Then the Green's function on each segment can be found in closed form as 

G^-(x,y,t) = 2-^[G^-(x,y,w)] 



0(t). 



Goc(x-y,t) 



F(x^y,t,ra/R) 



Gfn+(x,y,t) = 2 ^[G^+(x,y,a))] 



(51) 



(52) 



(53) 
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Gca{x + y,t) 



F{x + y,t,ra/R) 



G„±{x,y,t) = £ ^[G„±(x,y,ft))] 

^ 1 r 

C l2Rc} 



F{x + y,t,ra/R) 



where Goo(-^, t) is the Green's function of the passive infinite dendritic cable, 



Goo{x,t) = Z- 



g-y(ft))|jr| 

2Dy{w) 



1 



-t/r.-x^KADt) 



(9(f). 



(54) 



(55) 



(56) 



If an external stimulus Iia],m-ix, t) = 8{x — y)I{t) has a form of a rectangular pulse 
with I{t) = r\Q@{t)&{TR — t), the voltage response on each segment can also be 
found in closed form: 



V^-{xJ) = [B{x - y,t) - B{x - yj - tr) 

- {P(x + y,t)- P(x +y,t- tr))]/C, 
V„+(x, t) = [B{x +y,t)- B(x +y,t- tr) 

- (P{x + y,t)- P(x + y,t- tr))]/C, 
V„±(x, t) = [P(x + y,t)- P(x +y,t- tr)]/C, 



(57) 

(58) 
(59) 



where 



B(x,t) : 



^0 



-el^l/^erfc 



e-l^l/^erfcf- 



2JDt 



\2^Dt 



(60) 



P{X,t) : 



2DRgj 



[aF(x, t, ra/Rc]) + bF(x, t, ^/6) + cF(x, t, -^/s)], (61) 



and 



1 



Dt 



b = 



2v^(V^-r^//?Gj)' 



(ra/RGjV-S 



c = 



1 



(62) 



2v^(Ve + r^//?Gj)' 



These solutions generaUse earher results of Harris and Timofeeva [23] applicable 
to a neural network, but with gap-junctional coupling at tip-to-tip contacts of two 
branches. 



Appendix B: Proof of the "Sum-over- Trips" Rules at the Gap Junction 

As the Green's function is constructed in the Laplace domain, the rules at the gap- 
junction need to satisfy the following boundary conditions (after a length of each 
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Fig. 17 The GJ node with 
possible trips in its proximity 




branch, say labelled by k, of the tree is re-scaled sls X = yk(co)x, x e [0, Ck])'- 

Gm-j(0, Y, CO) = G^+j(0, Y, CO), 
Gn-j(0, Y, CO) = Gn+j(0, Y, CO), 



(63) 
(64) 



and 



+ 



dGjn+j 



X=0 



dx 



X=0 



= gGj{Gfn-j(0, Y, co) - Gn-j(0,co)), 

Ynjco) / dG^-j 
ra,n V 



dGn+j 



X=0 



■ gGj{Gn-j(0, Y, CO) - G^-j(0, Y, CO)). 



(65) 



(66) 



We prove here that the rules for generating the trip coefficients are consistent with 
these boundary conditions. 

Let X denote the distance away from the GJ node along the segment m~ (see 
Fig. 17). The location of the stimulus Y = yj(co)y, the segment number j and the 
variable co are all considered to be arbitrary. Suppose that we sum all the trips 
starting from the GJ node itself and ending at point Y on branch j. We denote 
the result of summing over all trips that initially leave the GJ node along segment 
m~ by Qm-j(^^ <^). along segment m+ by Q/^+j(0, Y, co), along segment n~ by 
Qn-jiO, Y, co) and along segment by Qn+j(0, Y, co). 

Trips that start out from X and move away from the GJ node are identical to trips 
that start out from the GJ node itself along segment m~ . The only difference is that 
the trips in the first case are shorter by the length X. We denote the sum of such 
shortened trips by Q^-j(—X, Y, co). The argument —X means that a distance X has 
to be subtracted from the length of each trip summed to compute Q (and not that the 
trips start at the point —X). 

Trips that start out from X by moving toward the GJ node and then reflecting 
back along segment m~ are also identical to trips that start out from the GJ node 
along segment m~ except that these are longer by the length X. In addition, be- 
cause of the reflection from the GJ node these trips pick up a factor —pG],n(o^) 
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according to the "sum-over-trips" rules. Therefore, the contribution to the solution 
Gm-j(X,Y,co) from those trips is —pG],n((^)Qm-j(^^^^(^)- Trips that start out 
from X by moving toward the GJ node and then continue moving along branch 
m, i.e. on segment m+, pick up a factor 1 — PG],n((^) and the sum of such trips is 
given by (1 — pG],n(co))Qm+ Finally, trips that start from X, move toward 
the GJ node and then leave the GJ node by moving out along segment n~ or 
pick up a factor PGJ,n((^) and contribute to the solution G^-j(X, Y, co) by the terms 

PGln(co)Qn-j(X, Y, co) OY PGln(co)Qn+ ji^^ 

The full solution G^-j (X, Y, co) includes the contributions from all different types 
of trips we have been discussing. Thus, 

G^-j(X, Y,co) = \ [Q Y,co) + {-pGln(co))Qm-j(X, Y,co) 

+ (1 - PGln(0)))Qm+j(X, Y, CO) + PGln(o))Qn- j (X , Y, Co) 
+ PGln(c^)Qn+j(X,Y,C0)l (67) 

The functions Q in this formula consist of infinite sums over trips, but we do not need 
to know what they are to show that the solution G^-j(X, Y, co) satisfies the GJ node 
boundary conditions. At the GJ node, we have 



Gm-j(0, Y,co) = 



1 



-[{l-pGln(co))Qm-j(0, Y, Co) 



DjYjico)' 

+ (l-m,^M)Qm+;(0, Y, CO) 

+ PGin(c^){Qn-j(0, Y, CO) + Q^+,(0, Y, co))]. (68) 

Considering the solution G^+j (0, Y, co) instead gives us the same expression as in 
(68) and, therefore, Gf^-j(X, Y, co) obeys the boundary condition (63). Similarly, we 
can show that 



G^-j(0, Y, co) = G^+j(0, Y, CO) 
1 



[{l-pGlm(co))Qn-j(0, Y, Co) 



DjYjico)' 
+ (l-m,m(^))Qn+;(0, Y, CO) 
+ PGlm(co){Q^-j(0, Y, CO) + Q^+j(0, Y, co))l 

which satisfies the boundary condition (64). 

To prove the boundary condition (65) we use Eq. (67) to find that 



(69) 



dG 



m J 



dx 



1 



dQm-j(-X, Y, CO) 



- PgjA^) 



dx 

dQm-j(X, Y,C0) 



dx 



^ springer 
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(1 - PGlnico)) 



dQm+j(X, Y, CO) 



dx 

dQn-j(X, Y, w) 



+ 



dQn+jiX, Y, w) 



dx 



dx 



X=0/ J 



x=o 



Using the following properties for the term Qkj (X, Y, co), k e {m ,m^,n ,n~^[ 
dQkjiX, Y, w) 



dX 

dQkj(-X, Y, w) 



dx 



= -Qkj(X, Y, CO), 
= Qkj(X, Y, CO), 



(70) 

(71) 
(72) 



Eq. (70) can be simplified as 

^ X=0 



dx 



(1 + PGln(c^))Qm-(0, Y, C0)-{1- /7GJ,n (^)) Qm+ (0, Y, Co) 

- PGin(co){Qn-(0, Y,co) + Q^+(0, Y,co)). (73) 



Similarly, 



dx 



X=0 



(1 + m,n(^))Qm+(0, Y, C0)-{1- PGln(co))Qm-(0, Y, Co) 

- PGin(co){Qn-(0, Y, co) + Q^+(0, Y, co)) . (74) 



Substituting (73) and (74) together with (68) and (69) in Eq. (65) gives us the right 
equality. Similarly, we can prove the boundary condition (66). 
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